In [1]:
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
mpl.rcParams['legend.fontsize'] = 10
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import leastsq
from math import acos, atan2, cos, sin
from numpy import array, float64, zeros
from numpy.linalg import norm
In [2]:
#%matplotlib inline
In [31]:
I = 10 #amps - change this back to 1.5
mu = 4*np.pi*1e-7 #This gives B in units of Tesla
In [19]:
R = .026 #meters
length = 0.25 #meters
guess_c1 = -1.26841572
guess_c2 = -5.81781983
guess_c3 = 0.04515335
p = np.linspace(0, 2 * np.pi, 5000)
z = p*length/(2*np.pi)
dp = p[1] - p[0]
In [20]:
def get_theta(c_1, c_2, c_3):
return c_1*p + c_2*p**2 + c_3*p**3
In [21]:
def get_x(c_1, c_2, c_3):
return R * np.cos(get_theta(c_1, c_2, c_3))
In [22]:
def get_y(c_1, c_2, c_3):
return R * np.sin(get_theta(c_1, c_2, c_3))
In [23]:
def get_z():
return p*length/(2*np.pi)
In [24]:
def cart2pol(x, y):
rho = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)
return(rho, phi)
In [25]:
def cart2pol(x, y):
rho = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)
return(rho, phi)
In [26]:
def j(g):
"""Returns numbers for list for line"""
l = 20.0*g + .2 #change y-intercept back to .2
return l
In [27]:
def B(zprime, c_1, c_2, c_3):
"""Returns B field in Tesla at point zprime on the z-axis"""
ex = get_x(c_1, c_2, c_3)
why = get_y(c_1, c_2, c_3)
r = np.vstack((ex,why,z-zprime)).transpose()
r_mag = np.sqrt(r[:,0]**2 + r[:,1]**2 + r[:,2]**2)
r_mag = np.vstack((r_mag,r_mag,r_mag)).transpose()
dr = r[1:,:] - r[:-1,:]
drdp = dr/dp
crossterm = np.cross(drdp,r[:-1,:])
return abs(mu*I/(4.0*np.pi) * np.nansum(crossterm / r_mag[:-1,:]**3 * dp,axis=0))
In [28]:
#this function was added to condense the return of function B into just a magnitude from its components
def Bmag(zpoints1, c_1, c_2, c_3):
Bdata = [1e4*B(zpoint, c_1, c_2, c_3) for zpoint in zpoints1]
Bdata1 = np.asarray(Bdata)
Bdata2 = []
for i in range(0, 5000):
Bdata2.append(np.sqrt(Bdata1[i,0]**2 + Bdata1[i, 1]**2 + Bdata1[i, 2]**2))
Bdata3 = np.asarray(Bdata2)
return Bdata3
In [30]:
zpoints = np.arange(0,0.15,0.00003)
k = []
for i in zpoints:
k.append(j(i))
d = np.asarray(k)
plt.plot(zpoints,d)
optimize_func = lambda points, c: Bmag(points, c[0], c[1], c[2])
ErrorFunc = lambda c,points,dat: dat[1000:3600] - optimize_func(points, c)[1000:3600]
c_initial = (guess_c1, guess_c2, guess_c3)
est_c, success = leastsq(ErrorFunc, c_initial[:], args=(zpoints,d))
print est_c
c1 = est_c[0]
c2 = est_c[1]
c3 = est_c[2]
Bdata = Bmag(zpoints, c1, c2, c3)
plt.plot(zpoints,Bdata)
ax = plt.gca()
ax.axvspan(0.03,0.12,alpha=0.2,color="green")
plt.ylabel("B-field (G)")
plt.xlabel("z (m)")
plt.show()
In [ ]:
Bdata = Bmag(zpoints, c1, c2, c3)
plt.plot(zpoints,Bdata)
ax = plt.gca()
ax.axvspan(0.03,0.12,alpha=0.2,color="green")
plt.ylabel("B-field (G)")
plt.xlabel("z (m)")
plt.plot(zpoints,d)
plt.show()
In [17]:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(get_x(c1, c2, c3), get_y(c1, c2, c3), z, label='solenoid')
ax.legend()
ax.set_aspect('equal')
plt.show()
In [ ]:
rho, phi = cart2pol (get_x(c1, c2, c3), get_y(c1, c2, c3))
plt.plot(-z, phi, '.')
plt.show()
In [56]:
#Why is this cell necessary?
plt.plot(get_theta(c1, c2, c3))
Out[56]:
In [9]:
#Calculate r vector:
r = np.vstack((x,y,z)).transpose()
plt.plot(r)
Out[9]:
In [10]:
# Calculate dr vector:
dr = r[1:,:] - r[:-1,:]
plt.plot(dr)
Out[10]:
In [11]:
# Calculate dp vector:
dp = p[1:] - p[:-1]
plt.plot(dp)
Out[11]:
In [12]:
# or the smart way since p is linear:
dp = p[1] - p[0]
dp
Out[12]:
In [13]:
r_mag = np.sqrt(r[:,0]**2 + r[:,1]**2 + r[:,2]**2)
plt.plot(r_mag)
Out[13]:
Converted the for loops to numpy array-based operations. Usually this just means taking two shifted arrays and subtracting them (for the delta quantities). But we also do some stacking to make the arrays easier to handle. For example, we stack x y and z into the r array. Note, this uses dp, and x,y,z as defined above, all other quantities are calculated in the loop because r is always relative to the point of interest.
In [14]:
def B2(zprime):
r = np.vstack((x,y,z-zprime)).transpose()
r_mag = np.sqrt(r[:,0]**2 + r[:,1]**2 + r[:,2]**2)
r_mag = np.vstack((r_mag,r_mag,r_mag)).transpose()
dr = r[1:,:] - r[:-1,:]
drdp = dr/dp
crossterm = np.cross(drdp,r[:-1,:])
return mu*I/(4*np.pi) * np.nansum(crossterm / r_mag[:-1,:]**3 * dp,axis=0)
In [15]:
B2list = []
for i in np.arange(0,0.15,0.001):
B2list.append(B2(i))
In [16]:
plt.plot(B2list)
Out[16]:
In [17]:
def B(zprime):
B = 0
for i in range(len(x)-1):
dx = x[i+1] - x[i]
dy = y[i+1] - y[i]
dz = z[i+1] - z[i]
dp = p[i+1] - p[i]
drdp = [dx/dp, dy/dp, dz/dp]
r = [x[i],y[i],z[i]-zprime]
r_mag = np.sqrt(x[i]**2 + y[i]**2 + (z[i]-zprime)**2)
B += mu*I/(4*np.pi) * np.cross(drdp,r) / r_mag**3 * dp
return B
In [18]:
Blist = []
for i in np.arange(0,0.15,0.001):
Blist.append(B(i))
In [19]:
plt.plot(Blist)
Out[19]:
In [20]:
Blist_arr = np.asarray(Blist)
B2list_arr = np.asarray(B2list)
In [21]:
plt.plot(Blist_arr - B2list_arr)
Out[21]:
In [ ]: